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Abstract 

We numerically explore the pasta structures and properties of low-density nuclear matter without any assumption on 
the geometry. We observe conventional pasta structures, while a mixture of the pasta structures appears as a metastable 
state at some transient densities. We also discuss the lattice structure of droplets. 
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1. Introduction 

In low-density nuclear matter, which is realized in the 
supernovae core or in the crust of neutron stars, inhomo- 
geneous structures of nuclear matter are expected JH. 
With increase of density, which ranges from well below 
to the normal nuclear density, the shape of constituent 
nuclei is expected to change from spherical droplet to 
cylindrical rod, slab, cylindrical tube, spherical bubble, 
and to uniform. These shapes are figuratively called as 
"nuclear pasta". 

The density-dependence of the species and the sizes 
of pasta are determined by minimizing the total energy 
density, i.e. the sum of the bulk, the surface, and the 
Coulomb energy densities. We can roughly assume that 
nuclear matter at sub-saturation density consists of a 
dilute gas phase and a dense liquid phase in chemical 
equilibrium, which determines particle densities in both 
phases. Thus the bulk energy density is independent of 
the structure, once the volume fraction is given, but the 
shape and the size of the structure are determined by the 
balance between the surface and the Coulomb energy 
densities. 

Among the early studies of nuclear pasta, geometri- 
cal symmetry of the structure was very often assumed, 



'Con'esponding author 

Email addresses: okamotoSnucl . ph . tsukuba .ac.jp 
(Minoru OKAMOTO), maruyama . toshikiSj aea . go . jp (Toshiki 
MARUYAMA), yabanaOnucl . ph . tsukuba .ac.jp (Kazuhiro 
YABANA), tatsumiOruby . scphys . kyoto-u .ac.jp (Toshitaka 
TATSUMI) 



i.e. the Wigner-Seitz (WS) cell approximation was em- 
ployed: The whole space is divided into equivalent cells 
with charge neutrality, and a geometrical symmetry is 
assumed with a given dimensionality. Then the shapes 
of the cell becomes sphere in three dimension (3D), 
cylinder in 2D, and plate in ID cases. Reflective bound- 
ary condition is imposed to the density distribution of 
each particle. There is no interaction between cells due 
to the charge neutrality. So all the physical quantities of 
matter are represented by those of a single WS cell. Fur- 
thermore, the calculation is only in one-dimension due 
to the symmetry, which drastically reduces the compu- 
tation. 

Within this approximation and assuming uniform 
density distributions in each phase, analytic expressions 
of the surface and the Coulomb energies can be derived 
for a given dimensionality of the symmetry. In 1983 
Ravenhall et al. have shown J2[] that the configuration to 
give the lowest energy changes from spherical nuclear 
droplet (3N) to cylindrical nuclear rod (2N), nuclear 
slab (1NB), cylindrical bubble (2B), spherical bubble 
(3B) and to uniform (U) with increase of the density. In 
1984 Hashimoto et al. have included the Coulomb inter- 
action among cells and have got essentially the same 
results. In this work the Coulomb energy was evaluated 
without the WS approximation, but only simple shapes 
of nuclei were assumed. 

On the other hand, Williams and Koonin have calcu- 
lated nuclear pasta using the Thomas-Fermi model for a 
system in a cubic cell with periodic boundary conditions 
In this calculation no assumption was made for the 
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structure of matter and they got essentially the same 
results as the studies with the WS cell approximation. 
There is also a recent calculation with the Hartree-Fock 
theory in a cubic cell with periodic boundary conditions, 
giving almost the same results again U- Though the 
above calculations did not assume any particular struc- 
ture, the size of the cell was rather small so as to include 
only one period of the pasta structures. The usage of 
small cell brings about an implicit but strong restriction 
to the matter structure, suppressing the appearance of 
complex structures. 

In this article we perform three-dimensional calcula- 
tions in periodic cubic cells with sufficiently large sizes. 
We will discuss how the pasta structures appear in the 
ground state of matter, showing their crystalline struc- 
tures. We also show some metastable states of matter 
which may be realized at finite temperatures. 

2. Method 

To describe the baryon interaction, we employ the 
relativistic mean-field (RMF) model with the Thomas- 
Fermi approximation |6]. The RMF model deals with 
fields of mesons and baryons introduced in a Lorentz- 
invariant way. It is rather simple for numerical calcula- 
tions, but realistic enough to reproduce the bulk proper- 
ties of nuclear matter. From the variational principle we 
get the coupled equations of motion for the mean-fields 
and the Coulomb potential as 

- V 2 cr(r) + »£o-(r) = go- N (p s p (r) + p s n (r)) 

-bm N gl N cr(rf + cg^crir) 3 (1) 
-V 2 w (r) + ml <o (r) = g<oN(p P (r) +p„(r)) (2) 
-V 2 R (r) + m 2 p R (r) = g pN (p P {r) - p„(r)) (3) 

V 2 V Co ui(r) = e 2 (p p (r)+p e (r)) (4) 

where pj(r) = (i^,(r)i/f,(r)), i = p,n is the nucleon scalar 
density. We use the same parameters as in Ref. ((J so 
as to compare the EOS and structural changes of pasta 
with or without WS cell approximation. 

Equations of motion for fermions simply yield the 
standard relations between the densities and chemical 
potentials, 

Pn = y]k F , n (r) 2 + m* N (r) 2 

+gojN<^o(r) - gpNRoir) (5) 
P-p = y]k F , p (r) 2 + m* N (r) 2 

+gcaN^o(r) + gpNRo(r) - Vcoui(r) (6) 
pAr) = -(Me-V Co ui(r)) 3 /3n 2 (7) 



where m* N (r) = - go-Nc(r) is an effective mass. We 
assume here the global charge neutrality. In case of cold 
catalyzed matter, the relation p„ - p p + p e is further 
imposed. 

To numerically simulate infinite matter, we use a cu- 
bic cell with periodic boundary conditions. We divide 
the cell into three-dimensional grid points. The best cell 
size should be as large as to include several periods of 
the pasta structures, and the best grid width should be as 
small as to attain smooth meson and fermion fields. Due 
to the limitations in the memory and the CPU time, we 
use the cell size of ~ 60 fm and the grid width \dr\ ~ 0.8 
fm. If only one or two periods of structure appear in the 
calculation cell, its shape may be affected by the bound- 
ary condition and the appearance of some structures is 
implicitly suppressed. 

Giving average densities of baryons (ps) and elec- 
trons, initial density distributions are randomly pre- 
pared. Then proper density distributions and meson 
fields are searched for. To obtain the density distribu- 
tions of baryons and electrons we introduce the local 
chemical potentials p„(r) {a = p, n, e). The equilibrium 
state is determined so that the chemical potentials are 
independent of the position. An exception is the region 
with an empty particle density, where the chemical po- 
tential of that particle becomes higher. We repeat the 
following procedures to attain uniformity of the chem- 
ical potentials. A chemical potential yu,(r) of a baryon 
i — p, n on a grid point r is compared with those on the 
six neighboring grids r' -r + dr, (dr - ±dx, +dy, ±dz). 
If the chemical potential of the point under considera- 
tion is larger than that of another pi(r) > /U,(r'), some 
part of the density will be transferred to the other grid 
point. This adjustment of the density is done simulta- 
neously on all the grid points. In addition to the above 
process, we adjust the particle densities between distant 
grid points chosen randomly so as to avoid making sep- 
arate droplets with different yu,. 

The meson fields and the Coulomb potential are ob- 
tained by solving Eqs. ((TJ-© using the baryon density 
distributions p,(r) (i = p, n) and the charge density dis- 
tribution p p (r) + p c (r). These equations are solved nu- 
merically by a conjugate gradient method. The electron 
density p c (r) is directly calculated from the Coulomb 
potential Vcoui(r) and the electron chemical potential p e . 
The global charge neutrality is then attained by adjust- 
ing fi e . Above processes are repeated many times until 
we get convergence. We used PRIMERGY BX900 of 
JAEA massively parallel computing system. To finish 
a calculation of a typical case, about 10 CPU days is 
needed. 
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3. Result 

We present here our first results with fixed proton 
fraction Y p for Y p - 0.1,0.3 and 0.5, leaving catalyzed 
matter in another paper. Shown in Fig. Q]are the proton 
density distributions in cold symmetric matter (proton 
fraction Y p = 0.5). We can see that the typical pasta 
phases with rod, slab, tube, and bubble, in addition to 
the spherical nuclei (droplet), are reproduced by our cal- 
culation in which no assumption on the structures was 
used. 




Figure 1: (color online) Proton density distributions of the ground 
states of symmetric matter (Y p = 0.5). Typical pasta phases are ob- 
served: (a) Spherical droplets with a fee crystalline structure at baryon 
density pg = 0.01 far 3 , (b) Cylindrical rods with a honeycomb crys- 
talline structure at 0.024 fnT 3 . (c) Slabs at 0.05 for 3 , (d) Cylindrical 
tubes with a honeycomb crystalline structure at 0.08 fin -3 . (e) Spher- 
ical bubbles with a fee crystalline structure at 0.09 fm~ 3 . 
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Figure 2: Density profiles of proton (red), neutron (green) and electron 
(blue) for ¥,, = 0.5. 

Any intermediate structures does not appear as a 
ground state at any density. In a droplet, we have seen 
that the proton density is highest near the surface due to 
the Coulomb repulsion, while the neutron density dis- 
tribution in the droplet is flat. Note that baryon density 
outside the droplets is zero for Y p = 0.3 and 0.5. The 
electron density is finite over all space but slightly lo- 
calized around the droplets. We can see this behavior 
of fermions for Y p = 0.5 and ps = 0.016ftrT 3 in Fig. 
[2] where plotted are the densities of proton, neutron and 
electron along a line which crosses through the droplets. 



We show the density dependence of the energy, the 
total pressure, and the baryon partial pressure in Fig. [3] 
The density dependence of these pressures is qualita- 
tively the same as the one with the WS cell approxima- 
tion. The difference between our results and those with 
the WS cell approximation, in which the same RMF 
framework is used, is the density region of each pasta 
structure; density region of the rod is wider and the tube 
narrower in our calculation. Figure [4] shows the radius 
of droplets Rj, the lattice constant R ce \\, and the volume 
fraction u of droplets. Here, Rj and R ce \\ are defined as 
follows, 
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and u = (Rd/Rce\\) 3 , where V denotes the cell volume, 
Nj the number of droplets in the cell, and the bracket 
(...) the average over the cell volume. Comparing the 
present results and those with the WS cell approxima- 
tion, the volume fraction shows the same behavior, but 
the lattice constant and the radius of the droplets are dif- 
ferent. 

When a spherical nucleus becomes large, it becomes 
unstable against fission. It is expressed as the Bohr- 
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> 2E sm f, where E 
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the Coulomb energy of a nucleus. In the WS cell ap- 
proximation, the Coulomb energy in a cell is expressed 
using the Coulomb energy of a nucleus as Ecoui ~ 
^Coui _ 3m^ 3 /2). On the other hand, the condition of 
optimum nuclear size gives E sm f - 2£coui- From these 
equations, the appearance of non-spherical nucleus in 
nuclear matter due to the fission instability has been ex- 
pected for the volume fraction u > 1/8. However, in our 
calculation, the structural change from droplet to rod oc- 
curs around a volume fraction u = 0.2 (see the value at 
Pb ~ 0.02 fm -3 in Fig. 0J. The relation between the 
Coulomb energies of a cell and that of a nucleus has 
been derived by using a uniform background electrons 
and uniform baryon density in a nucleus. The effect of 
the screening by charged particles, which is naturally 
included in both the WS cell calculation in Ref. Ja] and 
our present calculation, may be one of the origins of this 
difference. Also the difference of the droplet surface 
may be the origin, since they used the compressible liq- 
uid drop model with a sharp surface, while our droplets 
have a diffuse surface. 

There are some differences between our calculation 
and the one with the WS cell approximation. Calcu- 
lating in a large cell which has several periods of pasta 
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Figure 3: (color online) From the left, energy, pressure, and baryon partial pressure of symmetric matter. Red lines indicate droplet, green rod, blue 
slab, magenta tube, cyan bubble, and solid black uniform, respectively. Dotted lines are the EOS in the case of single phase 



structures, the interactions among the unit structures are 
properly included. By comparing both results, we can 
see almost the same behavior of the volume fraction, but 
slightly larger in our case. 
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Figure 4: Density dependence of the radius, the lattice constant and 
the volume fraction. Red lines are our results and green lines are those 
with the WS cell approximation 



One point unlike the conventional results emerged in 
the crystalline structure of droplets. In our calculation, 
it emerges as a face-centered cubic (fee) lattice, while it 
has been regarded to take a body-centered cubic (bec) 
lattice in the previous studies [8j,|9J]. Crystalline struc- 
tures in the bec and fee lattices give rise to a subtle dif- 
ference of the Coulomb energy, which amounts to about 
0.2-0.8 MeV. The ratio of the Coulomb energy for total 
energy difference is about 20%. In other words, most 
of energy differences come from the bulk energy of the 
nuclear matter. 
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Table 1 : Density dependence of droplet radii for the fee and bec lat- 
tices. 



Table[TJshows the baryon-density dependence of radii 
of droplets in the cases of the fee and bec lattices. The 
radii of droplets are different even if their baryon den- 
sities are the same. In Refs. it is reported that 
the bec crystalline structure of droplets is realized in 
the ground state at low densities, while the fee crys- 
talline structure appears in our calculation. This dif- 
ference may partially come from that they compare the 
bec and fee crystalline structures using droplets of the 
same size, while in our case they have different droplet 
sizes. In the QMD calculations ifiol [Till that precede 
the present calculations without assuming structure not 
only for static but also for dynamical aspects of nuclear 
matter, droplets form the bec crystalline structure. This 
difference might come from the treatment of electrons: 
In the QMD calculation, uniform electron distribution 
has been assumed, while in our calculation, inhomoge- 
neous electron distribution is attained consistently. We 
have performed another calculation with uniformly dis- 
tributed electrons. However, the result did not change 
and the fee structure is more favored. Thus it is con- 
firmed that the charge screening by electrons does not 
very much affect the crystalline structure |6fl. This dif- 
ference of the crystalline structure between the QMD 
calculation and the present calculation remains as a fu- 
ture problem. 

In Figs. and [6] we show the density dependence of 
the energy, the total pressure and the baryon partial pres- 
sure for Y p = 0.3 and 0.1, which are roughly relevant to 
the supernova matter and the neutron star crust. We ob- 
tain the typical pasta structure as ground states for any 
proton fraction above 0.1. Also in the cases of Y p -0.3 
and 0.1, the fee lattice of droplets is energetically more 
favorable than the bec. For rod and tube, the simple lat- 
tice is more favorable than the honeycomb lattice. In 
Fig. [7] we plot the density profile of proton and neutron 
for Y p = 0.1 with baryon density 0.016ftrT 3 along a line 
which crosses through the droplets. The neutron den- 
sity is finite at any point: the space is filled by dripped 
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Figure 5: (color online) Same as Fig.f3]for asymmetric matter with Y p = 0.3. 
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Figure 6: (color online) Same as Fig.f3]for asymmetric matter with Y p = 0.1. 



neutrons, which is in contrast to the case with higher Y p 
where the neutron density is zero outside the nucleus. 
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Figure 7: Density profiles of proton (red) and neutron (green) and 
electron (blue) for Y p =0.1. To show clearly the distribution of elec- 
tron, the electron density is magnified 20 times. 




Figure 8: Proton density distributions with complex structures (Y p = 
0.5). (a) mixture of droplet and rod, 0.022 far 3 , (b) slab and tube, 
0.068 fm" 3 ; (c) dumbbell like structure, 0.018 fm" 3 ; (d) diamond 
like structure, 0.048 fm" 3 . 



On the way of searching for the ground-state struc- 
tures, we sometimes observe exotic structures which are 
energetically metastable. Figure [8] (a) shows a mixed 
structure of droplet and rod at ps = 0.022 fm" 3 . Simi- 
larly, (b) is a mixed structure of slab and tube at 0.068 
fm" 3 . These structures appear around densities where 



the structures change. If these structures had appeared 
as ground states, transition from droplet to rod, and from 
slab to tube would happen more smoothly by way of 
those intermediate structures. In fact, the QMD cal- 
culations have reported some intermediate structures of 
droplet and rod, slab and tube as ground states 11 Oil . 

We have observed more exotic structures: Panel (c) 
at 0.018 fm" 3 looks like dumbbells which have been re- 
ported to appear in the dynamical compression of matter 
1 121] . Finally, panel (d) at 0.044 fm" 3 is a diamond struc- 
ture, which resembles double-diamond structure studied 
with a compressible liquid drop model J9[[T3|. These 
complex structures are difficult to be the ground states 
since they have larger surface area than simple pasta 
structures. At finite temperatures, however, these struc- 
tures might contribute to the Boltzmann ensemble. 

4. Summary 

We have numerically explored inhomogeneous struc- 
tures and properties of low-density nuclear matter using 
the RMF model and the Thomas-Fermi approximation. 
Without any assumption on the geometry, we have car- 
ried out fully three-dimensional calculations on large 
cubic cells with periodic boundary conditions. With 
increase of density, which ranges from well below to 
the normal nuclear density, we have observed that the 
ground state of matter shows the typical pasta struc- 
tures. More complex structures like "diamond" struc- 
ture, "dumbbell" structure, and mixtures of two types 
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of pasta structures appear as metastable states at some 
transient densities. As for the crystalline structures for 
the droplet and the bubble structures, fee lattice has been 
more suitable than bec lattice, which is different from 
the previous studies. For neutron star crust and super- 
novae matter, we should explore these structures and 
properties for/J-equilibrium and extend to finite temper- 
ature. At finite temperatures, complex structures might 
contribute to the Boltzmann ensemble. We can apply 
these contributions to the mechanism of glitch and cool- 
ing process of neutron stars and thermal and mechanical 
properties of supernovae. To explore jS-equilibrium mat- 
ter, we need more calculation space. But there are some 
problems in calculation performance. As mentioned be- 
fore, to finish a calculation of a typical case, many CPU 
time is needed. Most of the CPU time is consumed by 
the part to get the uniformity of chemical potentials and 
its iteration. To reduce the CPU time, we must improve 
the method and calculation code of this part. More real- 
istic calculation can be done by, for example, including 
gradient terms of the density to improve the description 
of the surface properties 11411 . 
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